function NE = NE_compute(positionk,fNE)
%UNTITLED3 此处提供此函数的摘要
phi=positionk(1);
lam=positionk(2);
h=positionk(3);
phi0=deg2rad(30.5279788680);
lam0=deg2rad(114.3557804579);


%%
R=rmnGet(phi);
RM=R(1);
RN=R(2);
N=(phi-phi0)*(RM+h);
E=(lam-lam0)*(RN+h)*cos(phi);
NE=[N;E];
fprintf(fNE,'%14.10f %14.10f\n',NE(1),NE(2));
end

%%
function R = rmnGet(phi)
%获得某个纬度处的子午圈半径：RM、卯酉圈半径：RN;
    a=6378137.0;
    e=0.08181919104;
    RM=(a * (1 - e^2)) / (1 - e^2 * sin(phi)^2)^(3/2);
    RN=a / sqrt(1 - e^2 * sin(phi)^2);
    R=[RM; RN];
end
